## This program is public domain
## Author: Paul Kienzle
function pdf = mvnpdf (x, mu = 0, sigma = 1)
  ## Check input
  if (!ismatrix (x))
    error ("mvnpdf: first input must be a matrix");
  endif

  if (!isvector (mu) && !isscalar (mu))
    error ("mvnpdf: second input must be a real scalar or vector");
  endif

  if (!ismatrix (sigma) || !issquare (sigma))
    error ("mvnpdf: third input must be a square matrix");
  endif

  [ps, ps] = size (sigma);
  [d, p] = size (x);
  if (p != ps)
    error ("mvnpdf: dimensions of data and covariance matrix does not match");
  endif

  if (numel (mu) != p && numel (mu) != 1)
    error ("mvnpdf: dimensions of data does not match dimensions of mean value");
  endif

  mu = mu (:).';
  if (all (size (mu) == [1, p]))
    mu = repmat (mu, [d, 1]);
  endif

  if (nargin < 3)
    pdf = (2*pi)^(-p/2) * exp (-sumsq (x-mu, 2)/2);
  else
%[vec,val]=eig(sigma);
%val(val<0)=eps;
%sigma=vec*val*vec';
    r = chol (sigma);
    pdf = (2*pi)^(-p/2) * exp (-sumsq ((x-mu)/r, 2)/2) / prod (diag (r));
  endif
endfunction
